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A  Parallel  Divide  and  Conquer 
Algorithm  for  the  Generalized  Real 
Symmetric  Definite  Tridiagonal 
Eigenproblem 

Carlos  F.  Borges*  and  William  B.  Gragg * 


Abstract. 

We  develop  a  parallel  divide  and  conquer  algorithm,  by  extension,  for  the  generalized 
real  symmetric  definite  tridiagonal  eigenproblem.  The  algorithm  employs  techniques 
first  proposed  by  Gu  and  Eisenstat  to  prevent  loss  of  orthogonality  in  the  computed 
eigenvectors  for  the  modification  algorithm.  We  examine  numerical  stability  and  adapt 
the  insightful  error  analysis  of  Gu  and  Eisenstat  to  the  arrow  case.  The  algorithm 
incorporates  an  elegant  zero  finder  with  global  monotone  cubic  convergence  that  has 
performed  well  in  numerical  experiments.  A  complete  set  of  tested  matlab  routines 
implementing  the  algorithm  is  available  on  request  from  the  authors. 


1  Introduction 

We  consider  the  problem  of  finding  a  matrix  U  €  5ftnxn  such  that 


Ut(T-SX)U  =  A-IX, 

is  diagonal,  or  equivalently 


UTSU  =  /  and  UtTU  =  A, 


0) 


where 


'Authors  address:  Code  Ma/Bc,  Naval  Postgraduate  School,  Monterey,  CA  93943.  Email: 
borges4fwayIon.math.nps.navy.inil 

•Authors  address:  Code  Ma/Gr,  Naval  Postgraduate  School,  Monterey,  CA  93943.  Email: 
gragg4tguinness.inath.nps.navy.nul 


1 


2  Carlos  F.  Borges  and  William  B.  Gragg 


’  01 

'  $1  7i 

ft  <*2  ft 

7i  ft  72 

ft 

,  s  = 

72 

ft-1 

7n  — 1 

.  ft»-l  ®n 

.  7n-l  ft 

and  S  is  assumed  to  be  positive  definite.  This  generalized  eigenvalue  problem  haa 
two  special  cases  that  are  of  interest  in  themselves.  They  are: 

1.  5  =  /,  the  ordinary  tridiagonal  eigenproblem. 

2 .  S  =  I  and  o*  =  0,  the  bidiagonal  singular  value  problem  (bsvp),  by  perfect 
shuffle  of  the  Jordan  matrix 


0 

B 


Bt 

0 


with  B  upper  bidiagonai  [16]. 

There  are  two  phases  to  the  divide  and  conquer  algorithm,  the  divide  (or  split) 
phase,  and  the  conquer  (or  consolidate)  phase.  We  shall  address  these  in  order. 


2  The  algorithm 

2.1  The  divide  phase 

Denote  by  ei  the  ith  axis  vector  where  the  dimension  will  be  clear  from  the  context. 
Let  s,  1  <  s  <  n,  be  an  integer,  the  split  index,  and  consider  the  following  block 
forms: 


7i 

e«-ift-i 

<>» 

ftef 

eift 

ft 

e,-i7«-i 

7»-ier_i 

ft 

7.ef 

ei7* 

53 

Note  that  s  =  n  is  possible;  then  T2,  S2,  and  ej  are  empty  [9,  10].  Suppose  we 
solve  the  subproblems 


U?  (7*  —  5*A)  £4  =  A*  -  IX  (k  =  1,2). 


(2) 
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2.2  The  conquer  phase 

The  conquer  phase  consists  of  solving  the  subproblems  (2)  from  the  divide  phase, 
consolidating  the  solutions,  and  finally,  solving  the  consolidated  problem.  Let 

ui  =  Uj  e,_i,  U2~l/?elt 

where  the  Ut  are  solutions  to  (2).  Then 

Ai  -  /A  -7,_iA) 

Ut(T-SX)U=  (A_,-7.-jAK  a, -6, A  (0,  -  7,A)u£  . 

u9(/?.-7.A)  A2-/A 

The  right  side  is  the  sum  of  a  diagonal  and  a  Swiss  cross: 


This  can  be  permuted  to  an  arrow  matrix  by  a  permutation  similarity  transfer- 
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M*"'  Aij-  H7"1'  *»]■ 

Since  S  and  j  uTq  ^,U  j  are  congruent  the  latter  inherits  positive  definiteness 
from  the  former.  Its  Cholesky  decomposition  is 

with  p2  —  7  —  uTC2u  >  0  the  Schur  complement  in  S  of 

[*  *]• 


d-i  _  I  -Cu/p 

1  i*  J 

and  a  second  congruence  transformation  with  R~ 1  gives 


A(X)  :=  R~TA(X)R~l 

=  r~t  f  ?„ 


D 

fiu 

u  7  B 

a 

w 

-  /A 

w 

R~x  -  IX 


(B-  DC)  u 

w  =  - , 

P 

o,  -  uT  (2 B  -  DC)  Cu 

W  =  - n - • 

P 2 

We  have  reduced  the  conquer  step  to  the  problem  of  solving  an  ordinary  eigen- 
problem  for  a  symmetric  arrow  matrix.  If  V  is  an  orthogonal  matrix  with 


AV  =  VA 


and  A  diagonal,  then  (1)  holds  with 


U  =  UP,R~1V 

'  Ui  -Uiu^.i/p  ' 

=  I/P  V. 

U2  -U2u2y,/p 
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It  is  useful  that  v*  =  t4u*  can  be  computed  in  O(n)  time  by  solving  SjVj  = 
e,_i  and  S2V2  =  ei  using  the  Cholesky  factorization  Si  =  L\Lj  and  the  reverse 
Cholesky  factorization  S2  =  L^L-i-  In  the  case  that  only  the  eigenvalues  are 
wanted  it  is  only  necessary  to  compute  the  first  and  last  rows  of  the  ^/-matrices 
which  constitutes  a  further  savings. 

In  summary,  the  conquer  phase  proceeds  by  consolidating  the  subproblems  and 
building  a  full  eigenproblem  for  an  arrow  matrix. 


3  Solving  the  eigenproblem  for  the  arrow 

In  this  section  we  consider  the  solution  of  the  eigenproblem  for  a  real  symmetric 
arrow  matrix 


A  = 


D 

bT 


b  ' 
1  . 


where  A  €  9?nx"  is  symmetric,  D  =  diag(a),  a  =  [01,  ...,o„_i]T,  ai  >  a2  >  ...  > 
a„_j,  and  b  =  [/3j ,  ...,/3„_i]t  >  0.  When  A  arises  from  the  bsvd  then  a  is  odd 
and  b  is  even,  that  is  a-f  7a  =  0  and  b  =  J b,  with  J  the  counter-identity,  the 
identity  matrix  with  its  columns  reversed,  aud  7  =  0. 

If  any  /3j  —  0  then  it  is  possible  to  set  A y  =  a;  and  deflate  the  matrix  since 
ej  is  clearly  an  eigenvector  (28).  We  shall  call  this  /3-deflation  and  note  that  if 
fij  <  to/^||b||  where  lolp  is  a  small  tolerance  then  a  numerical  deflation  occurs. 
We  derive  a  precise  value  for  tola  in  section  4.4. 

A  second  type  of  deflation  occurs  if  applying  a  2  x  2  rotation  similarity  trans¬ 
formation  in  the  (j,  j  +  l)-plane  that  takes  /3;  to  zero  introduces  a  sufficiently  small 
element  in  the  (j,  j+I)  position  of  the  matrix.  This  will  be  called  a  comio-deflation 
(see  [15]).  At  each  consolidation  step  we  perform  a  sweep  to  check  for  /^deflations 
followed  by  a  sweep  to  check  for  comiodeflations.  The  com fto-deflation  can  be 
arranged  so  that  the  ordering  of  the  a  j  is  preserved  whenever  one  occurs.  After 

deflation  the  new  /3*+i  :=  +  /?;2+,  >  /?J+ j  and  hence  no  further  /3-deflation 

can  occur.  The  co»n6o-deflations  can  be  disposed  of  with  a  single  pass  by  backing 
up  a  single  element  whenever  one  occurs.  Note  that  deflation  is  backward  stable 
since  it  results  in  small  backward  errors  in  A.  Deflation  for  the  BSVD  is  more 
delicate  involving  a  simultaneous  sweep  from  both  ends  of  the  matrix.  Care  must 
be  exercised  at  the  center  of  the  matrix. 

After  deflation  the  resulting  matrix  can  be  taken  to  have  all  0j  >  0  and  the 
elements  of  the  arrow  shaft  distinct  and  ordered,  that  is  <*1  >  a2  >  ...  >  q„_). 
An  arrow  matrix  of  this  form  will  be  called  ordered  and  reduced.  Henceforth,  we 
shall  assume  A  is  of  this  form. 

The  block  Gauss  factorization  of  A  -  XI  is 
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.4  -  A/  = 


1  0  I  f  D~XI 

bT(D-\I)~l  lj[  0r 


where  /,  the  spectral  function  of  A,  is  given  by 


b 

-/(A)  . 


n-1 


/(A)  =  A-7+£ 
i=i 


This  is  a  rational  Pick  function  with  a  pole  at  infinity  [1],  The  most  general  form 
of  a  rational  Pick  function  is 


/(A)  =  a-7  +  £-^-T>  6>  0.  (3) 

In  relation  to  the  various  divide  and  conquer  schemes,  the  case  6  >  0  corresponds 
with  extension,  6  =  0  with  modification,  and  6  =  y  =  0  with  restriction  [7]. 

Inspection  of  the  graph  of  the  spectral  function  reveals  that  the  elements  of 
the  shaft  interlace  the  eigenvalues 

A]  >  ft]  >  A2  >  ...  >  on-i  >  (4) 

Moreover,  in  the  present  case,  the  derivative  of  the  spectral  function  is  bounded 
below  by  one  so  that  its  zeros  are,  in  a  certain  sense,  well  determined. 


3.1  The  zero  finder 

The  fundamental  problem  in  finding  the  eigenvalues  of  an  arrow  is  that  of  providing 
a  stable  and  efficient  method  for  finding  the  zeros  of  the  spectral  function.  We 
now  examine  this  problem  in  some  detail. 

The  zero  finding  algorithm  we  present  is  globally  convergent  in  the  sense  that 
the  iteration  will  converge  to  the  unique  zero  of  /  in  (»i,a*_i)  from  any  starting 
value  in  the  closed  interval  where  we  put  ao  =  +00  and  q„  =  —00. 

The  zero  finder  converges  monotonicaily  at  a  cubic  rate  and  applies  to  a  general 
Pick  function  as  given  in  formula  (3). 


3.2  Interior  intervals 

The  iterative  procedure  for  finding  the  unique  zero  of  /  in  one  of  the  interior 
intervals  (04,0*- 1)  proceeds  as  follows.  Let  xo,  njt  <  z 0  <  a*-!,  be  an  initial 
approximation  to  A*.  If  rj  is  known  choose 
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=/(,)(*i).  ‘=0,1,2. 


Thus  ffjWo,  and  wj  must  satisfy 


1 

(fti-l  -  Xj)-1 

(ft*  -  Xj)'1 

a 

r  /(*,-) 

0 

(ft*-!  -  Xj)"2 

(fti  -  Xj)'2 

Wo 

=  /'(*>) 

0 

(ft*-!  -  Xj)-3 

(ftfc  -  Xj)'3  _ 

_  Wi  J 

n*j) 

and  we  find 


n  /  .  ,  .  0i  O, -fti-ia, -art 

a  =  3x,  —  (7  +<u-i  +  t*k)  +  > - . 

1  '  ft*  -  X;  a;  —  X*  ft,  -  X. 

ijU-l.k  ‘  J  *  J  ‘  1 

=  K--*,)*  /1  +  £  ff,  , 

^  „rr,, («!-*/ 

=  s+(«L1«nf,+  ^ 

ftt-i  -  at  ^  -  *;>*  n-  -  x>  ; 

Since  u0  >  0  and  wj  >  0  it  follows  that  Oj  is  a  Pick  function.  Thus  <f>j  has  a 
unique  zeto  xJ+]  €  Uu-,  ftjt-i)-  Also 


-'0  >  /*i_j  >  0  ,  wj  >  $  >  0. 


The  error  function 


/<*)-«*)=»-<■,+»>+  r  +stm 

has  n  zeros,  counting  multiplicities.  There  are  »— 3  zeros  exterior  to  (ftt,  ftt_i)  and 
three  more  at  x; .  Thus  the  error  function  crosses  zero  exactly  once  in  the  interval 
(at,afc_i).  Hence  Xj+i  lies  between  Xj  and  Afc,  and  the  iteration  is  monotonically 
convergent  from  any  starting  guess  x0  €  [o/t,a,t-i]  as  claimed.  The  cubic  rate  of 
convergence  follows  from  (5). 

Successive  iterates  can  be  found  by  solving  quadratic  equations.  Rather  than 
solve  <t>j{x)  =  0  for  Xj  +  i  it  is  better  to  solve 

4>j(xj  -  A)  =  0 

for  the  increment  A  =  Xj  —  Xj+1.  Some  rearrangement  using  (5)  reduces  this  to 


«A2  +  /iA  -  /  =  0, 


with 
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* 


(c*k-\  -  *;•)(*>  -  at)’ 

0  =  — - —  +  /(*>)-  (8) 

\ru_i-Z;  Qk-XjJ 

When  shifts  of  the  origin  to  the  nearest  pole  [15]  are  used  then  one  of  a*_i  or  a* 
is  zero.  The  computation  of  0  =  0(xj)  should  account  for  the  fact  that  it  has  only 
simple  poles  at  and  a*. 

If  we  start  at  the  midpoint  of  the  interval,  z<j  =  (<**-!  +o*)/ 2,  then  we  always 
have  0  =  0(xj)  >  /'(x,)  >  1.  This  can  be  seen  by  noting  that  0(xo)  =  /'(z 0) 
and  that,  when  xu  >  At-  then  for  all  of  the  succeeding  iterates  f(ij)  >  0,  by 
monotonicity,  and  —  —  4-  is  negative.  If  x0  <  A*  a  similar  argument 

applies.  It  follows  that  the  increment  can  always  be  computed  stably  as 


A  = 


■2f/0 


1  + 


\A+W 


(9) 


« 


3.3  Exterior  intervals 

The  treatment  of  the  two  exterior  intervals  is  geovielrically  the  same  as  above. 
Again,  the  approximating  function  has  poles  at  the  endpoints  and  the  residues  at 
these  poles,  and  the  constant  term,  are  chosen  to  satisfy  (5).  We  present  the  case 
for  the  interval  (c*i,oo),  the  case  for  the  other  exterior  interval  being  similar.  Now 


with 


<>j(x)  =  u0x  -  <r  + 


Wi 


a]  —  x 


n  -  I  .y, 

,  ,  0<  °i~°t 

u-o  -  1  +  >  7 - T7 -  >  1. 

-  -  (?=*)*  *  *• 


The  inequalities  are  strict  unless  n  =  2.  Again  we  find  (6)  where  now 


«- 1 


a  =  — 


i  +  ETJ 


JL 


(x,-a,)-  r, 


0  =  fUj)  + 


Xj  -  t»l 

These  are  limiting  cases  of  (7)  and  (8)  (  introduce  another  pole  oo  >  arj  and  let 
ao  — •  +oo).  If  x0  >  Ai  then  /(A)  >  0  so  0  >  f  >  1  and  A  is  again  computed 
stably  using  (9).  We  obtain  global  monotone  cubic  convergence  as  before. 
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Contrary  to  the  algorithms  of  [11,  12,  15]  our  algorithm  is  well-defined  when 
starting  at  the  endpoints  of  the  intervals.  The  algorithm  of  [23]  can  start  at  the 
*  endpoints  but  has  only  quadratic,  convergence. 

To  guarantee  that  xo  >  Ai  we  take  xo  to  be  the  iterate  in  («i,+oo)  from  +oo. 
As  Xo  — »  +oc  the  approximate  Pick  function  tends  to 

*  a;  -  7  +  (10) 

Cr  j  —  x 

Our  actual  starting  guess  is  the  zero  of  (10)  in  (o^-foo): 

'  «i  +  2ru  +  \/(2^)*  +  ii1>ir-;  -  ^>«>- 

Xo  =  '  ,  Him-' 

ft  1  H - TtUili— mmmm m  -f  <  ft  J  . 

When  shifts  are  used  we  have  at  =  0. 


3.4  Orthogonality  of  the  eigenvectors 

It  is  essential  that  the  computed  eigenvectors  of  the  arrow  matrix  be  numerically 
orthogonal.  As  a  point  of  entry  into  the  further  analysis  of  the  algorithm  we  now 
examine  the  orthogonality  of  the  eigenvectors  following  [15]. 

Consider  the  divided  difference 


/(A.fi) 


/(A)  -  /(/')  _ 

A  -  /<  ~  (ftj  -  A)(r»;  -  n) 

1  +  bT(D- A/)-I(£>-/(/rlb. 


Note  that  /i  =  A  gives  /'( A)  =  1  +  jj(/3  —  A/)  1  b[|.*.  If  /(A)  =  0  then 


(11) 


* 


v(A)  = 


A  —  <*j 

'  (A/  -  D)~] b  ‘ 

i 

1 

is  an  eigenvector  of  the  arrow  it  itrix  A 
value  A,  and 


D  h 
l>r  7 


associated  with  the  eigen- 


u(A)  = 


<’(  A) 

yrw) 


is  the  normalized  eigenvector  whose  last  element  is  positive.  The  ordering  of  A 
implies  that  its  matrix  of  eigenvectors  can  be  taken  positive  below  and  on  the 
diagonal,  and  negative  above. 
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Let  /( Ao)  —  /(/.„)  —  0  with  An  /<0.  Thus  An  and  po  are  distinct  eigenvalues 
of  A.  The  eigenvectors  u(A0)  and  w(/i0)  are  orthonormal: 

«(Ao)7  u(po)  =  /( A0 ,  m)  ~  0. 

Let  A  and  ft  be  approximate  eigenvalues  in  the  sense  that 


A  -  At, 

oj  —  Au 


IM< 


6 

i  +  6’ 


l‘  -  i<u 

<• J  -  /'ll 


(12) 


Here  6  >  0  is  hopefully,  bui  nut  necessarily,  close  to  the  machine  unit  t.  Note  that 
(12)  is  equivalent  with 


1  +  6> , 


('j  ~  /' 

")  ~  fo 


1  +  6'. 


These  conditions  imply  that  the  approximate  eigenvectors  u(A)  and  u(p)  are 
nearly  orthogonal.  For  we  have 


V7W(/0u(  A)Tu(p) 


Since 


/(A,/0-/(Ao,/to) 

a- 1 

E 


7T!  (t‘;  -  -  ,,) 

n-  I 

E 


•h 


fr't  ("j  -  A)("j  -  i0 


(l  -  (Qj  ~  A)(oj  -  fi)  \ 

V*  (O;  -  Ao)(oj  -  Po)/ 

(6;  +  6'  +  6j6j)  . 


ft  +  <;+v;isT^  +  rrCp<2« 


then 


\Z7'(A)/'(7<)H(A)Tri(/r)  =  26hr(/J  -  A/)-‘0(D  - 

with  |0|  <  /.  Thus 


\Z/'(^)/'(/0!w(A)Tu(/i)|  =  26jj(D  —  A/)-,l»)]2|)(D  —  /,/)—  1  bj|2. 


and  so 


|«(A)tu(/i)|  <  26. 

Condition  (12)  is  stringent.  If  we  let  (3t  — .  0  then  it  is  easy  to  show  that  A 
can  have  an  eigenvalue  Au  =  An(,4)  =  ofc  +  U{iif.)\  (12)  then  requires  that  the 
approximate  eigenvalue  A  sat  isfies  a  bound  l 
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|A-  A0|  <«(*$), 

which  is  difficult  if  /?t/l!NI  is  only  somewhat  larger  limn  machine  precision,  say 
f3/4  Two  techniques  are  used  to  attempt  lo  satisfy  (12)  -  shifts  of  the  origin  [15], 
and  simulated  extended  precision  (sep)  arithmetic  [26,  14].  Condition  (12)  means 
that 

jA  -  Aol  <  <5rnin{A0  -  tu,<U-i  -  -M 
When  shifts  are  used  it  means  that  A  is  nearly  //( Ao) ■ 


4  Numerical  stability  of  the  algorithm 

We  now  give  a  partial  analysis  of  the  stability  of  this  approach  to  the  eigenproblem 
for  the  symmetric  arrow  matrix.  Observe  that 

fix)  =  _  ruo-^) 

The  following  inverse  eigenvalue  problem  [6]  is  important;  given  {tty}  and  {Ay} 
satisfying  (4),  find  {/ )j }  and  ->  so  that  A(.4)  =  {A;}.  This  problem  is  simply  solved 
by  computing  the  residues  of  the  partial  fraction  decomposition  of  /.  In  particular 


For  fixed  {ory},  the  elements  of  the  arrow  head.  { ;4 )  and  y,  are  explicitly  known 
functions  of  the  eigenvalues. 

Now  let  {Ay}  be  a  set  of  approximate  eigenvalues  of  A  satisfying  (4).  Then 


0l 

7 


_  rr=i(f|*-  - 
ny**(n*  ~  °j) 

n  n  —  I 

EV£»,. 

i=i  i- i 


W4  >  0), 


(13) 

(14) 


define  a  modified  matrix  A  with  A{,4)  =  {Ay ).  To  obtain  a  backward  error  analysis 
for  the  complete  eigenvalue  problem  we  bound  the  deferences  fit  —  fit  and  7  —  7. 


12  Carlos  F.  Borges  aiul  William  H  Gragg 


4.1  Error  analysis  for  the  Dongarra-Sorensen  condition 

We  give  an  error  analysis  using  I  in-  Dongarra-Sorensen  condition 

T  =  *’■*'  IM<*.  (**) 

“t  -  A j 

where  6  =  0(()  is  of  the  order  of  the  machine  unit,  simplifying  that  in  [6]. 
Rearrangement  of  (15)  gives 

Aj  -  nk  =  (Xj  -  or)(l  +  ij,*}. 

It  follows  that 


■H  =  lH  II<>+  Am)  =  :U  I  1  +  <5;,, 

J=\  \  J=1 


and 


with  the  k  and  b" t  at.  most  only  slightly  larger  than  the  6jtt-  Thus 


where  6"  =  (2(c)  is  only  slightly  larger  than  b. 
Now  (M)  becoi tier. 


■k  -  -4 

,k 


j- i 

witli  one  of  the  poles  of  /.  Thus 

n 

i*  -  it  < 

J  =  ! 

To  minimize  this  hound  we  chouse  nt^}  to  he  a  pole  of  /  closest  to  Xj.  Clearly, 
a*U)  =  aj  <•*•(»)  =  o„_|.  so 

It  -  ll  <  f>  -  “|)  +  |A,  -  «*(;,(  +  (a„-i  -  A„) 

For  1  <  j  <  u  a  closest  pole  to  Xj  is  either  o;  or  o;  _  i .  The  distance 


r 


S# 


IAJ  -  “Kill  =  <">"  {•*/  “  1  ~ 


( 
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* 


is  maximized  when  A j  is  the  midpoint  of  the  interval  {o;,Oj_i),  and  the  value  of 
the  maximum  is  (o;  +  «j_i)/'2.  Thus 


|7  -7l 


1  *  ^ 

<  6  (Ai  —  til)  -f  ^  y^(n;-l  **  f>;  )  +  (c'u-l  ^n) 


J  =  - 


<  £(At  —  A„)  < 


In  summary,  the  Dongarra- .Sorensen  condition  implies  small  relative  errors  in 
each  and  a  small  absolute  error  in  7.  Fur  tin-  hsvu  this  implies  small  element¬ 
wise  relative  errors  since  the  condition  7=7  =  0  is  enforced  by  A j  +  A„+J  _7  =  0 
(only  half  of  the  eigenvalues  are  actually  computed,  the  rest  follow  from  this  con¬ 
dition). 


4.2  Rounding  error  analysis  of  the  computation  of  /(A) 

The  choice  of  a  termination  criterion  depends  on  a  careful  rounding  error  anal¬ 
ysis  of  the  particular  manner  in  which  ttc  compute  /(A).  Let  {«,-},  {/?«},  and  7 
be  floating  point  numbers.  We  represent  A  as  the  ordered  pair  of  floating  point 
numbers  (cr, /<)  where  the  slnfl  a  is  a  pole  closest  to  A,  and  A  :=  cr  +  ft.  For  the 
exterior  intervals  we  have  it  =  oj  or  it  =  n„_  For  the  interior  intervals  it  can  be 
determined  by  evaluating  /  at  the  midpoint  and  checking  the  sign.  We  compute 
/(A)  as 


t.-i 


1: 


c — '  II,  —  II 


J=l  ) 

with  the  standard  operation  precedence  rule.-,  where 


o'  =  Oj  —  it  ami  7'  =  7  —  a. 

We  use  Wilkinson’s  notation:  fl(x  *  »/)  =  (./■  *  j/)(l  +  6)  with  |6 j  <  c/(l  +  c) 
and  c  =  2~*  the  machine  unit.  More  generally,  i  denotes  numbers  not  essentially 
larger  than  2" *  (27)  and  the  rounding  errors  t  satisfy  |£|  <  c. 

We  define 


//(i»j  -  A)  :=  //(o'  -  /1)  =  //((n;  -  (r)  -  /1). 

v  If  a  =  o*  then 


//(o(t  -  A)  =  -//  =  m-  -  A, 

with  no  rounding  error.  For  j  ^  it. 
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* 

and  since  o*  is  a  pole  closest  to  A  then  |  <  2.  Thus  all  terms  etj  —  A  are 

computed  with  small  relative  error--: 


-  A)  =  (n,  ~  A)(  1  +  .1A,),  |A,  j  <  t. 


(16) 


When  computing  /(A)  =  we  add  the  term  A  -  ->  =  (A  -  a)  -  (7  -  <r)  last. 

A  routine  error  analysis  using  (Iti)  and 

|A-‘)|<|/(A)|  +  ’£r-i 


r  =  i 


I";  ~  M 


to  eliminate  the  term  jA  -  *, |  from  the  error  hound  gives 

(n-l 

:«i/tA)|-Hff-A|+(«+5)V7 - 

fa  K  " 


A| 


which  implies 


(n-l 

k-  A|  +  (»!+.' >)£ 


M 


(17) 


4.,'J  Termination 

Our  goal  is  to  choose  a  termination  criterion  so  that  we  stop  when  A  is  as  close 
to  the  true  eigenvalue  Ac  as  po-v-ibl.  .  Let  /<  =  Ac  in  (11)  with  /(A*)  =  0.  Now 
a*  <  Ac  <  <u--i  Also  let  oc  <  A  <  1  hen  the  terms  o^  —  A  and  Oj  —  A* 

have  the  same  sign  and 


I A  —  Ad  < 


l/(A)| 


^  =  •  |  n  ;  -  A  |  j  u  ,  —  A  t 


(18) 


rr 


To  obtain  an  upper  bound  for  |A  -  A*|  we  need  an  upper  bound  for  |/(A)|  and  a 
lower  bound  for  the  denominator.  For  the  latter  we  have 


»l  —  I 


'  +  £ 


;Ti  I'1-'  ~  At|o,  ~  A*  I  ~  maxr  I°J  -  A*r 


(19) 


Let  us  determine  how  small  |/(A)|  is  when  A  is  the  rounded  representation  of 
A*.  This  is 
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and  we  have 


A  :=  cr  +  fl(m  )  =  rr  +/i*.(  1  +  6) 
=  At  -f  in-d  —  ,\k  +  (A*  —  a)/! 


|(T  —  XL  |  =  min  jiij  —  At  | 

j 


|/(A)|  =  |A-A,j  l  +  V - r-i - 

\  fTi  K  -h  ) 


=  Kcr— A,)A|  1  +  V  - - 

\  l'*>  -  A)K  -  Ai) 


II  —  l  j'J 


<t-a1.|  +  Y'— L 

In,  - 


From  (17), 


\W('m<<  l  Iff -At.|+|Ji- rrl +  (,,  +  «)  ir-i- 
\  JTt  K  -  Ai 

Since  At  -  fr  =  (A  -  u)/{  1  +  f>)  then 


l//(/(A))i<(  a|A-rr| +  ,„  +  ())  y  ~u^. 

\  ,TI  I'D  -  *i 

We  terminate  and  sol  At.  :=  A  when 


!//(/( -M)|  i  ( ‘<^|A  -  <rj  +  ( n  +  0)  y  - — -  j  . 

\  fr>  i"i  -  *1/ 

Inequality  (17)  also  holds  if  /(A)  and  //(/(A))  av.  mimhangod.  Thus 


|/(At)|<.  (0|At-ff| +  (:$„+  \7)y-JL- 

{  U  k  -  a*i 

From  (18)  and  (19) 


ok-  A; !  +  (:)»+  17) T*-)  r-^t 
|A*  -A|.|<,  max  |A*  -  „,| - ±-Jz‘  h/7i 
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Since  \a-  At)  <  |«r  —  Ai-|  +  |A,  -  At|  ami  |rr  —  A4.  |  <  maXj  |«y  —  A*  j  the  computed 
eigenvalues  satisfy 


< 


|A*  —  At |  <  (.'{»+ I7)<  max  |<ij  —  A*( 

<  (i{n  +  e)t||/l||2. 


(21) 


4.4  Error  analysis  for  the  Gu-Eisenstat  condition 

From  •)  —*)  =  , ( A;  —  A; )  ami  ( J 1 )  w<-  liml 

I*/  -  *1  <  M"  +  li)‘ll-'lb- 

We  have  noted  that  the  Uoiigaira- Sorensen  condition  (15)  is  stringent.  It  is 
natural  to  ask  for  small  «/>««/«/>  errors  in  tin-  /4 •  If  we  replace  6jt  by  Sjjt/Pk  in 
the  analysis  in  section  4.1  we  liml  that 


and 


1-4  -  -4l  -/'■  i<'"i  <  ^(1+0(0). 

are  implied  by  the  «.%/«/  <> .mlilimi 

-'4— ~7  =  fj.l: 

<1/  -  a, 

We  must  bound  6. 

From  (20) 


»-i 


n-J 


I  A*  -  At)  l  +  £ 


it\  h'j  ~  Aril »,  -  A,,) 


<  mi  |At.  -  tr[  -F  ^ 


4 


i  trtj  ~  A*  | 


with  m  =  3(»  +  C).  Using 


|A-.T|<|A-A,|  +  iA,  -<r| 


and  the  Gu-Eisenstat  inequality. 


I 

K  -  A| 


< 


_ 1 _ 

(<•;  ~  A)(o7  -  A*.) 


lAf  -  A| 

(<* j  “  A)(«y  -  At)’ 


v 


we  get 


t 
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|A*  -  A*|  <  »n  | A*  -  <t\  + 


r-iii  —  I 

Z-3-  1 


r;: 


n  —  i 
i 


(t* A*  !(<*,- A*  ) 


where  t  has  been  increased  to  </(I  —  m<  ). 
By  Cauchy’s  inequality. 


|A*  —  At-  {  < 


< 


ji  u 


I  At  -  <rj  + 


V  &:!  — S^n) 


1/2 


tm 


^|Ai-  -  tr | 


for  every  j.  The  arithmetic-geoimi  i  ic  m.  mi.  inequality  and  the  triangle  inequality 
yield 


Thus 


}At  —  At.  |  <  we 


< 


< 


llll 


-  ffi  +  I'D 


A* 


\M\ 
ib  J 


-  At  | 


!<**  ~  fl  0j  \ 
|At-CDl!|b||j 


■^'"'Ill'll 


If  mc||b||  <  tij  for  all  j,  then 


IV* I  <  dim  j|b||. 


and  consequently 


|*'4  -  <4|  <  (i n(n  +  0j<  ||b|j. 

Thus  lol d  is  tm.  If  t)t  <  3(u  +  0 }< |[/»t j  w<  replace  Jk  by  zero  and  accept  a*  as 
an  eigenvalue  witli  nortnalized  eigenvector 

The  computed  eigenvectors  of  A  ar.  i  d;,  n  to  In*  those  of  the  nearby  matrix  A. 
Because  of  (13)  and  (10)  they  are  ronqiiited  t < •  high  relative  precision  elementwise 
and  hence  are  numerically  orthogonal 
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